#Changed?
# Establish directories
# NB BRT Source Functions must be located in the working directory

in.dir = '/home1/99/jc152199/brt/data/'
setwd(in.dir)
out.dir = '/home1/99/jc152199/brt/data/'

# Load the gbm library to perform Boosted Regression Tree Analysis
necessary=c('gbm','SDMTools','random','sp','rgdal')
#check if library is installed
installed = necessary %in% installed.packages()
#if library is not installed, install it
if (length(necessary[!installed]) >=1) install.packages(necessary[!installed], dep = T)
#load the libraries
for (lib in necessary) library(lib,character.only=T)

# Load source code for BRT functions
source('brt.functions.R.cjsedit.r')

# Read in data to model

model.data = read.csv('/home1/99/jc152199/MicroclimateStatisticalDownscale/ToAnalyse/MicroMacroMinMaxASCII.csv',header=T)

# Build brt model

brt.gbm = gbm(formula = micro_max~AWAP_max+coastdist+fpcmean+fpcvar+roaddist+solar, distribution = "gaussian", data = model.data, n.trees = 3000, interaction.depth = 3, shrinkage = .01, bag.fraction = 0.5, cv.folds=10)

# Establish directory for ASCII files
	
ascii.dir = '/home1/99/jc152199/MicroclimateStatisticalDownscale/250mASCII/'

# Import ASCII's

fpcmean.asc = read.asc(paste(ascii.dir,'STATIC/','fpcmeanwtplusbuffer250m.asc',sep=''))
fpcvar.asc = read.asc(paste(ascii.dir,'STATIC/','fpcvarwtplusbuffer250m.asc',sep=''))
coasdtdist.asc = read.asc(paste(ascii.dir,'STATIC/','coastdistwtplusbuffer250m.asc',sep=''))
disttoroad.asc = read.asc(paste(ascii.dir,'STATIC/','disttoroadwtplusbuffer250m.asc',sep=''))
rad015.asc = read.asc(paste(ascii.dir,'SOLAR/rad015.asc',sep=''))
awapmax.asc = read.asc(paste(ascii.dir,'STATIC/AWAPmaxwtplusbuffer250m.asc',sep=''))

# Making spatial predictions from a BRT model
# Start by making a vector of the names of the ASCII grids
# Then rename the columns from model data to the same as the name of the ASCII file
# Then name the variables for the BRT model the same as both the model and ASCII grid names

grid.names = c('AWAPmaxwtplusbuffer250m.asc','coastdistwtplusbuffer250m.asc','fpcmeanwtplusbuffer250m.asc','fpcvarwtplusbuffer250m.asc','disttoroadwtplusbuffer250m.asc','rad015.asc')

names(model.data)[c(4:8,11)]=grid.names

variable.names <- c(names(model.data)[c(4:8,11)])

# Set the working directory the location where the ASCII grids are stored or the scan function won't work

setwd('/home1/99/jc152199/MicroclimateStatisticalDownscale/250mASCII/STATIC/')

for(i in 1:length(grid.names))
	
	{
	
	assign(variable.names[i],scan(grid.names[i], skip=6, na.string = "-9999"),pos=1)

	}

# Bind all the grid data into a single dataframe
	
predict.data = data.frame(AWAPmaxwtplusbuffer250m.asc,coastdistwtplusbuffer250m.asc,fpcmeanwtplusbuffer250m.asc,fpcvarwtplusbuffer250m.asc,disttoroadwtplusbuffer250m.asc,rad015.asc)
	
# Run the spatial prediction function

gbm.predict.grids(brt1, predict.data, want.grids = T, sp.name = "micro_max",pred.vec = rep(-9999,2095980), filepath = "/home1/99/jc152199/MicroclimateStatisticalDownscale/250mASCII/", num.col = 1086, num.row = 1930, xll = 254162.1, yll = 7812803, cell.size = 250, no.data = -9999, plot=T)

# Read in the spatial prediction

micro_max.asc = read.asc(paste(ascii.dir,'micro_max.asc',sep=''))

# Create a basic image of the spatial prediction

png(paste(in.dir,'micro_max.png',sep=""),width=600,height=1000,bg=rgb(.5,.5,.5,0,names='clear')) #Command bg=rgb(.5,.5,.5,0,names='clear') produces a transparent colour
      
par(pty="m",oma=c(5,0,5,0),mar=c(0,0,0,0),cex=1,xpd=T)
      
image(micro_max.asc,axes=F,ann=F,frame.plot=F,oldstyle =T,col=c(heat.colors(9)[9:1]))
      
dev.off()
	  







